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Abstract 

We apply an unsteady Reynolds-averaged Navier-Stokes (URANS) solver for the simulation 
of a synthetic jet created by a single diaphragm piezoelectric actuator in quiescent air. This config- 
uration was designated as Case 1 for the CFDVAL2004 workshop held at Williamsburg, Virginia, 
in March 2004. Time-averaged and instantaneous data for this case were obtained at NASA Lang- 
ley Research Center, using multiple measurement techniques. Computational results for this case 
using one-equation Spalart-Allmaras and two-equation Menter’s turbulence models are presented 
along with the experimental data. The effect of grid refinement, preconditioning and time-step 
variation are also examined in this paper. 


Introduction 

Significant interest has been growing in the aerospace community in the field of flow control 
in recent years. An entire AIAA conference is now devoted every other year to this field. In 
March 2004, NASA Langley Research Center, in conjunction with five other international orga- 
nizations held the CFDVAL2004 workshop, 1 in Williamsburg, Virginia. The primary objective of 
this workshop was to assess the state of the art for measuring and computing aerodynamic flows 
in the presence of synthetic jets. Thomas, Choudhari, and Joslin 2 have conducted an exhaustive 
and comprehensive survey identifying the feasibility of using active flow control to improve the 
performance of both external and internal flows. Suggested applications cover a wide range from 
smart materials and micro-electro-mechanical-systems (MEMS) to synthetic zero net mass jets for 
enhancing control forces, reducing drag, increasing lift, and enhancing mixing, to name a few. 
It is also conjectured that active flow control would permit the use of thicker wing sections in 
non-conventional configurations, such as the blended wing body (BWB) configuration, without 
compromising the aerodynamic performance. 
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Most of the research in the area of active flow control is of an empirical nature, mainly due 
to the cost and lack of confidence in computational methods for such complex flows. However, 
without the availability of efficient and well-calibrated computational tools, it will be a very diffi- 
cult, expensive, and slow process to determine the optimum layout and placement for active flow 
control devices in practical applications. With the continuous reduction of computer costs in re- 
cent years, researchers are devoting more attention to the simulation of such unsteady flows and 
flow control devices from a computational point of view [Ref. 3-8]. With few exceptions, most 
of the numerical studies are undertaken without an active interaction with experimental investiga- 
tors. Comparisons with experimental data are sometimes done years after the experimental data 
have been acquired. Under such a scenario, one has to reconstruct many of the details about the 
experimental arrangement and boundary conditions without the benefit of concrete and consistent 
information. Based on our experience from previous validation exercises, 9 we recognized the need 
for active collaboration of the computational and experimental research. Without a symbiotic re- 
lationship between the two groups, major misunderstandings can develop when results from these 
disciplines differ significantly. We were very fortunate to have a cooperative relationship with the 
researchers conducting the experiments, as well as access to pertinent experimental data. 

Our primary objective for this work is to calibrate an existing computational scheme with ex- 
perimental data for the time-dependent flows encountered in active flow control environments. We 
devote special attention to establishing appropriate boundary conditions for such flows, especially 
in the absence of the detailed experimental data required for closure. 

The configuration chosen for CFD validation is identified as Case 1 in the CFDVAL2004 work- 
shop 1 and represents an isolated synthetic jet formed by a single diaphragm, piezoelectric actuator 
exhausting into ambient quiescent air. Multiple measurement techniques, including PIV, LDV, and 
hotwire probes were used to generate a large body of experimental data for this configuration. Ref- 
erences 1 and 10 describe the details of the experimental setup and geometric configuration. In 
this paper, we assess the effects of grid refinement, preconditioning and turbulence models on the 
computational simulations of the flow field generated by this flow control device. We replace the 
actuator cavity with a simpler configuration in the present simulations. We demonstrate and cali- 
brate our computational method for simulating synthetic jets by comparing the numerical results 
with the experimental data. 


Governing Equations 

A generalized form of the thin layer Navier-Stokes (N-S) equations is used to model the flow. 
The equation set is obtained from the complete N-S equations by retaining the viscous diffusion 
terms normal to the solid surfaces in every coordinate direction. For a body-fitted coordinate 
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system (£, r/. () fixed in time, these equations can be written in the conservative form as: 


fl(U) (9(F — F v ) <9(G — G v ) i9(H - H v ) _ 

dt d £ dr] d( 


( 1 ) 


where U represents the conserved variable vector. The vectors F, G, H, and F v , G v , H v represent 
the convective and diffusive fluxes in the three transformed coordinate directions, respectively. 
In Eq. (1), Vol represents the cell volume or the Jacobian of the coordinate transformation. A 
multigrid-based, multiblock structured grid flow solver TLNS3D, developed at NASA Langley 
Research Center is used for the solution of the governing equations. References 1 1 and 12 describe 
the TLNS3D solver in detail, therefore only a brief summary of its general features is included here. 


Discretization 

The spatial terms in Eq. (1) are discretized using a cell-centered finite volume scheme. The 
convection terms are discretized using second-order central differences with a matrix artificial dis- 
sipation (second- and fourth- difference dissipation) added to suppress the odd-even decoupling 
and oscillations in the vicinity of shock waves and stagnation points. 13 ~ 15 The viscous terms are 
discretized with second-order accurate central difference formulas. 11 The zero-equation model 
of Baldwin-Lomax, 16 one-equation model of Spalart-Allmaras 17 and Menter’s two-equation SST 
model 18 are available in TLNS3D code for simulating turbulent flows. For the present computa- 
tions, the Spalart-Allmaras (SA) model and the Menter’s SST model are used. 

Regrouping the terms on the right-hand side into convective and diffusive terms, Eq. (1) can 
be rewritten as: 

^ = -C(U) + D p ( U) + D a ( U) (2) 

where C(U), D p { U), and D a (U) are the convection, physical diffusion, and artificial diffusion 
terms, respectively. These terms include the cell volume or the Jacobian of the coordinate trans- 
formation. 

The time-derivative term can be approximated to any desired order of accuracy by a Taylor 
series 

— = ^[aoU ra+1 «i un a 2U n 1 + asU n 2 + ...] (3) 

The superscript n represents the last time level at which the solution is known, and n + 1 refers 
to the next time level to which the solution will be advanced. Similarly, n — 1 refers to the solution 
at one time level before the current solution. Eq. (3) represents a generalized backward difference 
scheme (BDF) in time, where the order of accuracy is determined by the choice of coefficients 
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ao, a i, a 2 ... etc. For example, a 0 = 1.5, a\ = —2 and a 2 = .5, results in a second-order 
accurate scheme (BDF2) in time, which is the primary scheme chosen for this work because of 
its unconditional stability and good robustness properties. 19 Regrouping the time-dependent terms 
and the original steady-state operator leads to the equation: 


— — — T J n9_ i + ^ UV " 1 '") 

At At 


S{ U n ) 


( 4 ) 


where E(\J n,n ~ 1 ’ •■) depends only on the solution vector at time levels n and earlier. S represents 
the steady-state operator or the right-hand side of Eq. (2). By adding a pseudo-time term, we 
rewrite the above equation as: 


— + ^U” +1 + 
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Solution Algorithm 

The algorithm for solving unsteady flow relies on the steady-state algorithm in the TLNS3D 
code. 1112 The basic algorithm consists of a five-stage Runge-Kutta time-stepping scheme for ad- 
vancing the solution in pseudo-time. Efficiency of this algorithm is enhanced through the use of 
local time- stepping, residual smoothing and multigrid techniques developed for solving steady- 
state equations. Because the Mach number in much of the domain is very low, we consider the use 
of preconditioning methods 20 ’ 21 to improve the efficiency and accuracy of the flow solver. 

In order to solve the time-dependent Navier-Stokes equations (Eq. 5), we add another iteration 
loop in physical time outside the pseudo-time iteration loop. For fixed values of E(XJ n,n ~ 1 ’-), 
we iterate on U n+1 using the standard multigrid procedure of TLNS3D developed for steady- 
state, until the the pseudo-residual + 'ppp — ,S'(U) approaches zero. This strategy, originally 
proposed by Jameson 22 for Euler equations and adapted for the TLNS3D viscous flow solver by 
Melson et. al, 19 is popularly known as the dual time-stepping scheme for solving unsteady flows. 
The process is repeated until the desired number of physical time steps is completed. The details 
of the TLNS3D flow code for solving unsteady flows are available in Ref. 19,23,24 and 25. 


Boundary Conditions 

The boundary conditions required for solving the Navier-Stokes equations, such as the no-slip, 
no injection, fixed wall temperature or adiabatic wall, far-held and extrapolation conditions, are 
well understood and readily available in most how codes including the TLNS3D code. On the 
other hand, accurate simulation of oscillating diaphragm requires information about mode shapes 
and moving grid capability. The mode shape information is not readily available for this configu- 
ration, therefore we simulate this boundary condition by a periodic velocity transpiration condition 
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imposed at the diaphragm surface. The frequency of the transpiration velocity at the diaphragm 
surface in the numerical simulation corresponds to the frequency of the oscillating diaphragm. We 
obtained the peak velocity at the diaphragm surface from numerical iteration to match the exper- 
imentally measured peak velocity of the synthetic jet emanating from the slot exit. The pressure 
at the moving diaphragm is also required for closure. However, in the absence of unsteady pres- 
sure data from the experiment, we imposed a zero pressure gradient at the diaphragm boundary. 
We also tested the pressure gradient boundary condition obtained from a one-dimensional normal 
momentum equation. 26 This had very little impact on the solutions. Because of its simplicity and 
robustness, we selected the zero pressure gradient boundary condition at the diaphragm surface. 

Synthetic Jet Test Case: Background 

The test configuration is a single diaphragm piezoelectric actuator operating in quiescent air. 
The oscillatory motion of the diaphragm produces a synthetic jet that exhausts into the surrounding 
air. This configuration, shown in Fig. 1, consists of a 1.27 mm wide rectangular slot connected 
to a cavity with a piezoelectric diaphragm and corresponds to Case 1 of the CFDVAL2004 work- 
shop on flow control devices. 1 The cavity and diaphragm geometry of this actuator are highly 
three-dimensional in the interior. However, the actual slot through which the fluid emerges is a 
high aspect ratio rectangular slot and can be modeled as a two-dimensional configuration. Partial 
view of the 2-D grid provided to the workshop participants, and used by the present authors, is 
shown in Fig. 2(a) and 2(b). The computational results contributed by the workshop participants 
are available in Ref. 1 and 27. A consensus developed during the workshop that simulating the 
flow field inside the actuator cavity with an oscillating piezoelectric diaphragm from first principles 
was beyond the capability of existing CFD codes. Most of the workshop participants, including 
the present authors, modeled the internal cavity of the actuator as a two-dimensional configura- 
tion and simulated the diaphragm motion via a transpiration condition imposed at the diaphragm 
surface. Some of the workshop participants further simplified the cavity modeling by imposing a 
transpiration condition at the bottom part of the slot’s neck or even directly at the slot exit. After 
examining these results, we concluded that as long as the unsteady velocity signal at the slot exit 
replicates experimental conditions, details of the cavity modeling have an insignificant effect on 
the development of the synthetic jet emanating from the slot. Another conclusion derived from 
the CFDVAL2004 workshop was that no particular methodology or turbulence model emerged as 
superior for simulating this test case. 1 - 27 

We include here sample results obtained with the TLNS3D code to encapsulate the status of 
CFD simulations at the conclusion of the CFDVAL2004 workshop. For these computations, we 
used the 2-D grid of approximately 65,000 nodes from the CFDVAL2004 workshop website as 
the baseline grid. As depicted in Fig. 2(a) and 2(b), this grid includes the internal cavity and the 
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diaphragm geometry. A sinusoidal transpiration condition is imposed at the diaphragm surface, 
which is represented by the left face of the cavity in Fig. 2(b). We present the time history of 
the phase-averaged v-velocity at x = 0, y = 0.12 mm, and the time-averaged v-velocity along 
the jet centerline in Fig. 3 and 4, respectively. We performed the computations for three grid 
densities, two physical time steps and two turbulence models. The coarse grid (eg) was obtained 
by eliminating every other point from the baseline grid, and a fine grid (fg) was obtained by adding 
50% more points in the normal direction. We see from these figures that varying all these factors 
did not produce any major changes in the computational results in the region near the slot exit. 
However, further away from the slot exit, these results indicate that the coarse grid (eg) does not 
provide adequate resolution for this problem. A finer grid (fg) and a reduction in time step (low dt) 
had an insignificant effect on the computational results. However, significant differences are seen 
in outer region between the SA and the SST computations. It is difficult to make any definitive 
conclusion regarding the accuracy of the turbulence models based on these limited comparisons. 

One of the major difficulties identified during the CFDVAL2004 workshop was the large dis- 
parity in experimental data obtained from different measurement techniques, as seen in Fig. 3 
and 4. Such a variation in experimental data made it difficult to validate the numerical methods. 
Part of the difficulty in acquiring a consistent set of experimental data arose from the fact that the 
performance of the piezoelectric diaphragm depends on ambient conditions. Also, its performance 
degrades over time, which means that, for a given input voltage, the actuator produces smaller jet 
velocities as it ages. Because these experiments were conducted over several months, inconsisten- 
cies crept in the data. 


Results: Synthetic Jet II 

Yao et al. 10 have recently revisited the synthetic jet test case and acquired experimental data for 
this configuration with a new piezoelectric diaphragm. They obtained the detailed field data with 
the PIV technique and pointwise data along the jet centerline with hotwire and LDV techniques. 
They monitored the performance of the actuator regularly and demonstrated good consistency 
among multiple measurement techniques. 10 

Based on the CFDVAL2004 workshop 1 results, it can be concluded that replicating the flow 
conditions at the slot exit is more important than the detailed modeling of cavity geometry for 
accurate simulation of the growth of a synthetic jet for this configuration. Therefore, we simulated 
the new experimental test case with a simplified cavity geometry, shown schematically in Fig. 5. 
We imposed the transpiration condition at the bottom of the slot’s neck to simulate the velocity 
generated by the oscillating diaphragm. A similar boundary condition treatment has also been 
studied in detail by Yamaleev and Carpenter. 28 They demonstrated that for actuators with deep 
cavities, specifying the transpiration condition at a distance of at least 4-5 slot widths away from 
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the slot exit produces only a small loss in numerical accuracy. A top-hat velocity profile, with a 
dominant frequency of 450 Hz replicating the experimental conditions, was imposed at the bottom 
boundary. The precise form of the velocity signal was obtained by curve fitting the measured 
velocities at the slot exit (x = 0, y = 0.3 mm) with a fast Fourier transform (FFT) to reflect the proper 
mode shape and to ensure zero net mass transfer. The amplitude of this velocity was determined 
numerically to match the peak velocity from the experiment at the slot exit. The free stream Mach 
number in the exterior quiescent region is specified as M 0 0 = .001 to simulate incompressible flow 
in the compressible flow code to avoid numerical difficulties at Mach zero. 

The computational grid for this case was obtained from the baseline grid described in the 
previous section by eliminating the portions of the grid below the neck of the slot. The resulting 
grid consists of over 60,000 nodes and provides adequate resolution based on the grid refinement 
study reported in Ref. 1 and as seen from the results in Fig. 3 and 4. Similarly, 72 time-steps/period 
corresponding to a 5° phase angle between the time-steps provides sufficient temporal resolution. 
The baseline results were obtained with the one-equation Spalart-Allmaras (SA) turbulence model. 
In addition, the solutions on the baseline grid were also obtained with low-speed preconditioning 
(prec) and Menter’s turbulence model (sst). Finally, the results on fine grid (fg) with the SA 
model have been obtained. Based on the peak jet velocity and slot width, the Reynolds number 
is approximately 3000, and therefore falls in the regime where the jet is expected to be turbulent. 
Therefore, we assumed the flow to be fully turbulent in the present computations. 

The time history of the vertical velocity for a complete period from the computational results 
is compared with the experimental data in Fig. 6 at x = 0 and y = 0.3 mm. This is the closest point 
to the slot exit where the PIV data is available. In addition to the PIV data, LDV measurements are 
also available at this location and are shown in this figure. The LDV data was scaled down by a 
factor of 0.9 as suggested by Yao et al. 10 to match the diaphragm displacement for the two sets of 
measurements. Unlike the earlier results shown in Fig. 3, the overall agreement between the com- 
putational and experimental results is quite good at this location. The four sets of computational 
results shown in this figure are indistinguishable from one another, indicating a minimal effect of 
preconditioning, grid density and turbulence model at the slot exit boundary. 

In Fig. 7, we compare the time-averaged v-velocities along the jet centerline with the PIV 
and LDV data. The experimental data from two different techniques (PIV and LDV) are in fair 
agreement with each other. The overall agreement of the baseline TLNS3D results with the exper- 
imental data is also quite good. The low-speed preconditioning results included in this figure are 
essentially identical to the baseline results. Because low-speed preconditioning 20 - 21 ' 24 primarily 
reduces the artificial viscosity for unsteady flows, we may infer that the artificial viscosity is low 
in these simulations, even without preconditioning. Similarly, the effect of grid refinement (fg) 
is seen to be negligible on the computational results. However, the SST results differ from the 
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baseline results in regions away from the slot exit (y > 8 mm). 

Figures 8 and 9, respectively, show the time-averaged v-velocity profiles at y = 1 and 4 mm. 
Except for a smaller velocity peak at y = 1 mm, the computational results are in very good agree- 
ment with the experimental data and with each other at these locations. Comparing the contour 
plots of the time-averaged v-velocities obtained from measured PIV data and baseline TLNS3D 
computations (Fig. 10(a) and 10(b), respectively) gives a global perspective of the velocity field. 
Although the computational results are available over a much larger domain, these figures show 
a domain covering a distance of only 8 mm from slot exit, corresponding to the region for which 
the high resolution PIV data was available. The computational results accurately capture all of the 
prominent features seen in the PIV data including the width and spreading rate of the synthetic jet. 
The contour plots obtained with preconditioning and the SST model (not shown here) for the same 
domain are nearly identical to the baseline computations. 

We now examine the phase-averaged velocities at selected locations in space and time, starting 
with v-velocities at y = 2 and 4 mm along the jet center line. Figures 11(a) and 11(b) show the 
PIV and FDV data, along with TFNS3D solutions at these locations. The computational results 
are in fairly good agreement with the two sets of experimental data, especially in the suction 
phase. The agreement with the experimental data further away from the slot exit is slightly worse 
during the peak expulsion cycle. In particular, the CFD results predict a delayed phase shift for the 
peak expulsion, reflective of a smaller convective speed for outward movement of the synthetic jet 
compared to the experimental data. Except for a slightly larger peak velocity for the SST model 
during the expulsion phase, all four sets of computational results are practically indistinguishable 
from one another. 

We gain a broader perspective of the flow field by examining the contour plots of the velocities 
at the phase angles representative of the expulsion (phase = 75°) and suction (phase = 255°) cycles. 
Figures 12-19 show the velocity contours for streamwise (u-vel) and vertical velocities (v-vel) ob- 
tained from the PIV data and TFNS3D computational results (baseline). Although not shown here, 
the results obtained with preconditioning, the SST model and finer grid show very little variation 
from the baseline solutions over this domain. These figures were generated using identical contour 
levels for both the experimental and CFD data to provide quantitative comparisons. The solid lines 
represent positive values for the velocities while the dashed line represent negative values. This 
sign convention is helpful in identifying the flow direction and the position of the vortex center. It 
is clear from these figures that the computational results capture most of the pertinent features ob- 
served experimentally and are in very good agreement for the suction phase. During the expulsion 
phase, the computed vortex center is located closer to the slot exit compared to the experimental 
data, although the peak velocity at the vortex center is in good agreement with the PIV data. Yao 
et al. 10 have observed increasing three-dimensional effects for this case as one moves away from 
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the slot exit, mainly because of ring vortices formed from the slot edges. We conjecture that these 
ring vortices induce forces that accelerate the convection of the synthetic jet in the far field. 


Concluding Remarks 

Detailed comparisons have been presented for time-averaged and phase-averaged velocities 
between experimental data and CFD results. The effect of truncation errors were found to be 
small based on grid refinement, preconditioning, and physical time-step refinement studies. The 
development of the synthetic jet in the quiescent medium is driven primarily by the velocity field 
at the slot exit. Hence, approximating this forcing function is much more crucial than detailed 
modeling of the cavity and parametric variations of the numerical algorithm. The computational 
results in the reduced domain with a modified forcing function are found to be in much better 
agreement with the new experimental data in the near field. However, the agreement with the 
experimental data deteriorates in regions further away from the slot exit. Based on the available 
experimental data, it appears that the flow becomes three-dimensional after 5-6 slot widths away 
from the exit. Future work should focus on 3-D computations for this configuration to resolve such 
issues. 
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(a) Global view 


(b) Detailed view 


Figure 2 Computational grid for Piezoelectric Actuator 
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Exp. data, PIV 




Figure 3 Time-history of v-velocity at x = 0, Figure 4 Time-averaged v-velocity along jet 

y = 0.12 mm centerline 



Figure 5 Simplified model of actuator 
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Figure 6 Time-history of v-vel near slot exit Figure 7 Time-averaged v-vel along centerline, 

(x = 0, y = 0.3 mm), Synjet II Synjet II 



x (mm) 



x (mm) 


Figure 8 Time-averaged v-vel at y = 1 mm, Syn- Figure 9 Time-averaged v-vel at y = 4 mm, Syn- 
jet II jet II 
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V-vel (m/sec) 




(a) PIV measurements (b) TLNS3D computations 

Figure 10 Time-averaged v-velocity contours for Synjet II 




(a) x = 0, y = 2 mm 


(b) x = 0, y = 4 mm 


Figure 11 Phase-averaged v-velocity comparisons for Synjet II 
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Figure 16 Contour plots of u-vel, PIV data at 
phase = 255° 



x, mm 


Figure 18 Contour plots of v-vel, PIV data at 
phase = 255° 
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Figure 17 Contour plots of u-vel, TLNS3D re- 
sults at phase = 255° 
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Figure 19 Contour plots of v-vel, TLNS3D re- 
sults at phase = 255° 


17 OF 17 



